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Abstract 

We present an algorithm for doing Gibbs sampling on a quantum computer. The 
algorithm combines phase estimation for a Szegedy operator, and Grover's algorithm. 
For any e > 0, the algorithm will sample a probability distribution in 0{^=) steps with 
precision 0(e). Here 5 is the distance between the two largest eigenvalue magnitudes 
of the transition matrix of the Gibbs Markov chain used in the algorithm. It takes 
0{V) steps to achieve the same precision if one does Gibbs sampling on a classical 
computer. 
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1 Introduction 



In Ref. PQ, Szegedy proposed a quantum walk operator for each classical Markov 
chain. In Ref. [2], Somma et al. proposed a method for doing simulated annealing on 
a quantum computer. In Ref. [3], Wocjan et al. showed how to improve the Somma et 
al. algorithm. The algorithms of Somma et al. and Wocjan et al. both use Szegedy 
operators. In Ref. [3], I presented computer programs called QuSAnn and Multiplexor 
Expander that implement ideas of Refs.[2] and [3], and also some of my own ideas 
about quantum multiplexors. 

In Ref. [5] , I described one particular algorithm for doing Gibbs and Metropolis- 
Hastings sampling of a classical Bayesian network (i.e., a probability distribution) on 
a quantum computer. In this paper, I describe a different algorithm for doing Gibbs 
sampling on a quantum computer. Unlike my first algorithm, this one uses Szegedy 
operators. For any e > 0, this new algorithm will sample a Bayesian network in 
steps with precision 0(e). Here 5 is the distance between the two largest eigenvalue 
magnitudes of the transition matrix of the Gibbs Markov chain used in the algorithm. 
It takes O(^) steps to achieve the same precision if one does Gibbs sampling on a 
classical computer. 

This paper assumes that its reader has read the section entitled "Notation and 
Preliminaries" in Ref. [5]. The reader should refer to Refs.jHlH] for clarification when 
any notation of this paper eludes him. 

2 Dual Gibbs Markov Chains 

In this section, we will discuss dual "Gibbs" Markov chains with transition matrices 
Mi and M 2 , respectively. These two transition matrices are both defined in terms of 
a single classical Bayesian network x. 

2.1 Definitions of M 1 and M 2 

Consider a classical Bayesian net with N n d s nodes, labeled x l: x 2 , . . . ,x N where 
Xj G for each j. (As usual in my papers, I indicate random variables by underlining 
them.) Let x = (x l7 x 2 ,...,x N ). Let x assume values in a set Sx which has 
N s = 2 Nb elements. 
Let 

tt(x) = P(x = x) (1) 

for all x G Sx- 

For N n ds = 3 and x, y G Sx, let 

Mx{y\x) = Px 1 \x 2 ^(yi\x2,x 3 )P^\^ l (y 2 \x 3 ,yi)P^\^ 2 (y 3 \y 1 ,y 2 ) , (2) 
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and 



M 2 (y\x) = Px 1 \x 2 ^(yi\y2,y3)Px 2 \^x 1 (y2\y3,x 1 )P^ 1 ^ 2 (y 3 \x 1 ,x 2 ) . (3) 

(M 2 (y\x) can be obtained by swapping Xj and yi in the conditioned arguments of 
M 1 {y\x).) Note that J2 y M i(y\ x ) = 1 and Y. y M 2{y\x) = L Define M l and M 2 for 
arbitrary N nc i s using the same pattern. M\ and M 2 are transition matrices of the type 
typical for Gibbs sampling. (See Ref.[5] for an introduction to Gibbs sampling and 
the more general Metropolis-Hastings sampling). 

You can check that 7r() is not a detailed balance of either M\ nor M 2 separately. 
However, the following property is true. We will refer to this property by saying that 
7r() is a detailed balance of the pair (Mi, M 2 ). 

Claim 1 

M x {y\x)-K{x) = M 2 (x\y)n(y) (4) 

for all ijG Sx ■ 
proof: 

Let P(xj, Xk, ■ ■ •) stand for P(xj = Xj,x k = x k , ■ ■ ■)■ Assume N nds = 3 to begin 
with. One has 



Mi(y\x) P(yi\x 2 , x 3 )P(y 2 \x 3 , yi)P(y 3 |yi, 1/2) 
M 2 (x\y) P{xi\x 2 , x 3 )P(x 2 \x 3 ,y 1 )P(x 3 \y 1 , y 2 ) 
P(yi, x 2} x 3 )P(y 2} x 3 , yi)P(y 3 , y u y 2 ) 
P(x!, x 2 , x 3 )P(x 2 , x 3 , yi)P(x 3 ,yi, y 2 ) 
P(yi, x 2} x 3 )P(y u y 2l x 3 )P{y) 
P(x)P(y 1 , x 2 , x 3 )P{y x , y 2l x 3 ) 

P(y) 

P(x) ' 

A proof for an arbitrary number N n d s of nodes follows the same pattern. 
QED 

2.2 Eigenvalues of M u M 2 and M hyb 
Let 



(5) 
(6) 
(7) 



Ajiylx) = yjMj(y\x) , (9) 



for j = 1,2 and x, y G S^. It's convenient to define a hybrid function of Mi and M 2 , 
as follows: 
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M hyb (y\x) = A 2 (x\y)A 1 (y\x) (10) 

for x,y G S%- (Note that unlike Mi(y\x) and M 2 (y\x), Mh y b(y\x) is not a probability 
function in y, its first argument.) 
Define the quantum states 

|(7r)"} = 5>(a;)Hz} (11) 

a; 

for 77 = |, 1. (Note that only the r\ = | state is normalized in the sense of quantum 
mechanics.) 

Claim 2 

M J \n) = \n) for j = 1,2, (12) 

and 

M^ fe |v^) = |v^> • (13) 
Also, Mi, M 2 and M hyh have the same eigenvalues. 

proof: 

Taking the square root of both sides of the pair detailed balance statement 
Eq.(|4|), we get 



Ai(?/|x)v / ^) = A 2 {x\y)^{y) . (14) 

Therefore, 

M %6 (y|x) =A 2 (x|y)^i=A 2 (a;|y)v / ^) = -==M 2 (x\y) yffiy) . (15) 

Hence, 

£ M^x)^) = M 2 (x\y)ir(y) = ir(y) , (16) 

X X 

£ M 2 (a%)7r(y) = £ Mx^/lxMx) = vr(x) , (17) 

and 

J2M hyb (y\x)^/^ = J2^= M 2^\y)V^)V^= VW>- (is) 



Order the elements of the finite set Sx in some preferred way. Use this preferred 
order to represent Mi, M 2 and Mh y b as matrices. Define a diagonal matrix D whose 
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diagonal entries are the numbers ir(x) for each x G S%_, with the x ordered in the 
preferred order: 

D = diag[(-ir(x)\ x ] . (19) 

Since 

Ml = D~ X M X D , Ml yh = D-^M 2 Di , (20) 

it follows that 

(let (Mi - A) = det(M 2 - A) = det(M %6 - A) (21) 

for any A G C. 
QED 

Let the eigenvalue^] of Mh y b (and also of Mi and M 2 ) be m , mi, . . . mjv s -i G C 
with m = 1 > |mi| > jm^l ■ ■ ■ > |mjv s -i|- Define \rrij) to be the corresponding 
eigenvectors of M hyb (but not necessarily of Mi and M 2 ). Thus 

M hyh \m 3 ) = mjlrrij) , (22) 

for j — 0,1, ... , Ns — 1. In particular, \m ) = ly/n). 

For each j, define tpj G [0, |] and m G [0, 27r) so that rrij = e %r)i cos (fj. (Thus, 
cos (fj > 0). Note that m = 1 so ipo = 0. The Mi eigenvalue gap 5 is defined as 
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5 = 1 — |mi|. 5 ~ -y when 991 is small. 

3 Q-Embeddings C/i and 

In this section, we will define a "q-embedding" Uj of Mj, for j = 1,2. (For more 
information about q-embeddings, see Ref.[5].) 

For simplicity, we begin this section by considering a Bayesian net with only 
3 nodes x 1; a; 2 ,X3, and such that each of these nodes is binary (i.e., S^. = Bool for 
j = 1, 2, 3). At the end of this section, we will show how to remove these restrictions 
and make our treatment valid for general Bayesian networks. 

Using the same language as Ref . [5J , consider a unitary matrix U\ of the form 
shown in FigJTJ with its multiplexor gates defined as follows. Let x_Ak) G Bool 
and x[Ak) G Bool for any j,k. U\ has 3 analogous gates (a.k.a. nodes) labeled 
(x;(2),x 3 (2>,x 2 (2>), (x 2 (3),x' 1 (3),x 3 (3)), and (x' 3 (4), x 2 (4), a£(4». Consider the first 
of these for definiteness. Let the probability amplitude A(x[(2), £3(2), x 2 (2) ^(l), 23(1), ^2(1)) 
of node (x / 1 (2), ^3(2), x 2 (2)) satisfy the constraint 

1 There must be a single eigenvalue 1 and all others must have a magnitude strictly smaller than 
one because of the Frobenius-Perron Theorem. The eigenvalues may be complex. 
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Figure 1: Unitary matrix U\ expressed as a product of quantum multiplexors. 



A(x' 1 (2},x 3 (2),x 2 (2)\x' 1 (l) = 0,x 3 (l),x 2 (l)) (23) 
= ^J^^W(2)|x3(2) ) x 2 (2))C:; 2 1 > > C3 3 ( 2 1 )- (24) 
If we indicate non-zero entries by a plus sign, 
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(25) 




for some #g G R. Here the right pointing arrow means that the expression at the 
origin of the arrow can be extended to the expression at the target of the arrow. 
From the above definition of U\, it follows that, for x,x',y,y' G Bool 3 , 
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Hence, 
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Figure 2: Unitary matrix U2 expressed as a product of quantum multiplexors. 

Besides U\, it is convenient to consider another unitary matrix called L^- We 
define U2 to be of the form of Figj2j where the multiplexors are defined in such a way 
that U2 satisfies, for all x,x',y,y' G Bool 3 , 



(y\ rj |o> 
(y'\ U2 W) 
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(31) 
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Hence 



or U 2 \x')\0) m = |x')(A 2 |a; / )) . 



(32) 



\x') 



W) 



Uj is called the q-embedding of M, for j = 1,2. 



A = 




Figure 3: The unitary matrix A is a quantum embedding of a probability matrix 
P(b\a), where a has Nbo, bits and b has Ngb bits. 



So far we have considered the q-embeddings Ui and U 2 for the case of a classical 
Bayesian network x with 3 binary nodes. What if x has N n d s nodes and some of those 
nodes have more than 2 states? In that case, we must use several qubits (horizontal 
lines) for each node (and an equal number of qubits for the dual node a^) in 
FigsJH and [2j More specifically, suppose P(xi|x 2 , x 3 , . . . , xjq nds ) equals P{b\a) where 
a G Bool N}3 ^ and b G BooI Nb ^. For the number of bits Nsa, define the number of 
states N Sa = 2 Nb &. Likewise, let N Sb = 2 Nb *. The constraint Eq.flU]) generalizes to 



A(b,a\b = 0,a) = ^P(b\a)5* , (33) 

where a, a G BooI Nb ± and b, b G Bool NBb -. Eq.(|33]) can be expressed in matrix form 
as follows: 



[A(b,a\b = 0,a)] 



where, for all b G BooI Nb ^, D b '° G M. Ns&XNs& are diagonal matrices with entries 





(6 = 0,a) 


(6, a) 


D o,o 


1 


D i,o 







(34) 



(D b >°) ara = ^P~WK ■ 



(35) 
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By adding more columns to the matrix of Eq. (}34"|) . one can extended it (see section 
entitled "Q-Embeddings" in Ref.[5]) to a square matrix which can be expressed in 
terms of multiplexors as in Figj3j 

The Markov Blanket MB(i) for a node ^ of the classical Bayesian network x 
satisfies (see section entitled "Notation and Preliminaries" in Ref.JS]) 

P(xi\x {i} c) = P(xi\x M B(i)) ■ (36) 

If the set MB(i) is strictly smaller than the set {i} c , this property can be used to 
reduce the number of controls for the multiplexor in U\ and U2 corresponding to 

P(Xi\x {i} c). 

Given the two q-embeddings U\ and U2 for a Bayesian network x, we can define 
a unitary matrix U as follows 



U = UlUt . 

Matrix U has the following highly desirable property: 
Claim 3 For any j, k G {0, 1, . . . , N s - 1}, 



(m 



'3 v j 



proof: 



(0| jAjj \m k ) 

K-l u ' Ul 10) 



E 



y.x 





\ (y|Af ' 




\x) 


rrijly) 






Ai\x) 



(x\m k ) 



^(rnj\y)Kl{y\x)K 1 (y\x)(x\m k ) 
mj\M hyb \m k ) = rrijSj . 



y.x 



(37) 



(38) 



(39) 
(40) 
(41) 



QED 



4 Szegedy Quantum Walk Operator W 

In this section, we will define a special type of Szegedy quantum walk operator W 
corresponding to a Bayesian net x. We will then find the eigenvalues of W. 



4.1 Definition of W 

As in Ref.[l], define the projection operator n and its dual projection operator tt by 
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*- |0 - ' ioxoi ' < 42 > 

Then the Szegedy quantum walk operator W for the Bayesian net x is defined by 

W = U{-lfU\-lf . (43) 

4.2 Eigenvalues of W 

To find the eigenvalues of W, we will use the following identities. 
Claim 4 

7r|mj0) = \rrij0} , (44a) 
7t(U tylmjO) = rrijlmjO) , (44b) 
Tt(t f/t )l m i°) = rn*\ mj 0) , (44c) 

for allje{0,l,...,N s -l}. 
proof: 

From the definition of ft, we see that 

* !°» = | 0> , . (45) 



\rrij) \m,j) 



Also, 



MUt) l ° } =T |0)( °!c/! mi) =m- |0) (46) 
* \rrij) ' |0) J |mj) ' 



and 



♦tt^)! 0> » =e ! 0>< wn! ^ ! 0> % =™;!\ ■ w 

QED 

An immediate consequence of Claim H] is that 

(mjO\U X \m k 0) = ( mj 0\7tU | |m fc 0> = mfi) , (48) 

forj,fce{0,l,...,JV s -l}. 

Note that since m = 1, Eq. (|48|) implies that 



|m 0) = ?7 % |m 0) . (49) 
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Another consequence of Claim H] is that |moO) is a stationary state of W. 
Indeed, one has 



W\m 0) = U{-l)*U^{-lf\m 0) (50) 

= Ul(l-2n)lU\-l)\m 0) (51) 

= (l-2m £/£K-l)K0) (52) 

= (l-2)(-l)|m 0) (53) 

= |m 0). (54) 



Let 



^ s ^^n{| mj 0),[/|K0)} (55) 

for j G {0,1,..., Ns — 1}- (By "span" we mean all linear combinations of these 
vectors with complex coefficients.) 

Claim 5 WVl sy = V{ usy for all j e {0, 1, . . . , N s - 1}. For j = 0, let 

|-0o> = |m 0) . (56) 
(IV'o)} is an orthonormal basis for V busy and W\i])q) = \^o) . For j ^ ; let 

fe) = J* * |m ^' 0) ~ e± ^'K°» • ( 57 ) 

V2 sm yjj 

iV'-j')} ^ 5 an orthonormal basis for Vl usy and W\ip±j) = e ±t2<Pi \i})±j) . 
proof: 

Using the identities of Claim HI one finds after some algebra that 

W\mfi) = (-1)^-0) + (2m*)U I \rrij0) , (58a) 

and 

W(UX)\mfi) = (-2771^)1^^0) + (-1 + A\ mj \ 2 )UX \rrij0) (58b) 

for all j. 

According to Eqs. fl58|) . Vl usy is invariant under the action of W for each j. By 
virtue of Eq. (j4"8~]) . V busy is 1-dim for j = and 2-dim if j ^ 0. We've already proven 
that |raoO) is a stationary state of W. 

Now consider a fixed j ^ 0. Both U(—1) % W and (—l) n are reflections, and 
reflections are a special type of orthogonal matrix, so the product of these 2 orthog- 
onal matrices is also an orthogonal matrix. In fact, it's a rotation about the axis 
perpendicular to the planar subspace V 3 busy . The vectors \rnj0), and U \ \mfi)} are 
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|e 2j > A e iT W£lm ; 0> 



lm$=|e,.> 




Figure 4: Definition of \e±j) and \e2j). 



independent but not orthogonal. However, we can express them in terms of orthogo- 
nal vectors (see FigH]) as follows: 



\rrijO) = \e\j) 



(59a) 



and 



e~ ir »U t \rrijO) = cos(^0|ey) + sin(^)|e 2j ) . (59b) 
In the \e\j), \e2j) basis, we find after substituting rrij = e tVj cos (<£,■) into Eqs.f l58|) that 



W 



cos(2ifij) sin(2^) 
— sm(2<pj) cos(2ifij) 

The eigenvalues of this matrix are e <Pj , with corresponding eigenvectors: 



(60) 



(61) 



These eigenvectors satisfy 



(V-ilVi) = , (iPMj) = = 1 • (62) 

By expressing \e\j) and |e 2j ) in Eq.( l6"Tj) in the original basis, we get Eq. (157)1 . 
QED 

Define the following vector spaces: 



V = span{\x) <S> \y) ■ x,y G S^} 



(63) 



V A = span{\x) (g) |0) : x G S^} 



V B = U % V A 



(64) 
(65) 



and 



Vbusy = Va + Vb- 



(66) 
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V can be expressed as a direct sum of Vbusy and its orthogonal complement V^ : 

V = V busy © VL y ■ (67) 
From Claim El it follows that Vbusy is a direct sum of the subspaces V^, : 

iVs-l 

V** = VLy • (68) 
i=o 

Recall that matrices Mi , M2 and M/^ are Afg dimensional whereas W is iVj di- 
mensional. Since the size of S% is Ng, dim{V) = N$. From Eq. (1681) and Claim El 
dim(Vbusy) = 2N S - I. Furthermore, {{i/ij) : j = 0, ±1, ±2, . . . , ±(N S - 1)} is an 
orthonormal basis for Vbusy 

At this point we've explained the action of W on Vbusy, but we haven't said 
anything about the action of W on V^ SJ/ . Next we show that W acts simply as the 
identity on V^ sy . (This is what one would expect since the vectors in V^ sy are parallel 
to the axis of the W rotation.) Recall that if S and T are subspaces of a vector space 
V, then (S + T) L = S ± D T ± . Therefore, 

= vi n vi . (69) 
From the definitions of Va and Vb, it's easy to see that 



and 



Claim 6 



for all \4>) EV, 



Vi = span{\x) ® \y) : x G and y G - {0}} , (70) 

v£ = tf$(vi). (71) 

= |0) (72) 



proof: Let \<j>) G V^ sy = Vj[ H V#. Hence |0) G Vi and |0) = C/ | |0) for some 
\9) G Vi. 



l)*Z7 t (— 1^ 



t/t(-lft^(-l)°l0) 

u $ (1 - 2tt) $ tfu t \e) 

Ut(l-2n)\9) 



(73) 
(74) 
(75) 
(76) 



QED 
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It's interesting to compare the present paper with Ref.jl]. For Ref. jl], M\ = 
M 2 = M and tt() is a standard detailed balance for M instead of a detailed balance 
for the pair (M U M 2 ). For Ref.g], M hyb = M sym , m, = m*, XJ 1 = U, U 2 = U, 
U = U\U X = U^U, U =1 W X- When U =$ W X as in Ref.pEj, Eq.flEEj) and Eq.( li4c]) 
are essentially identical, whereas in the M\ ^ M 2 case, it's less obvious that these 
two equations are true simultaneously. 



5 Quantum Gibbs Sampling Algorithm 

In this section, we will describe an algorithm for doing Gibbs sampling on a quantum 
computer, utilizing the Szegedy operator W that we have so painstakingly discussed 
in previous sections. 

We begin by choosing^) some xo G Sx for which P(x = xq) ^ 0. Now define 

|x 0> = \x = x ) ® \0)® Nb . (77) 

Note that \xqQ) G Vb usy and 

(ipo\x$) = (Vtt\x = xo) = y/n(x ) . (78) 

y/ ir(x) = a/ P{x) can be easily evaluated at a single point x = x . Our quantum 
Gibbs algorithm consists of performing the original Grover algorithm with beginning 
state \xq0) and target state \ipo). Define the following 2 reflection operators 

R beg = (-I)l»o0><<o0| ^ (79) 

and 

Rtar = (-l)l^>^l . (80) 

RbegRtar is a rotation by an angle O(^tt(x )) in space span{ \i/jq), \x 0)} C Vbusy Let 

L = 0{-1==). (81) 

If a/ 7t(x ) = Oil/ yJNs), then L iterations of Rb eg R tar will take the beginning state 
to the target statejj To implement this use of Grover's algorithm, we need to compile 
(with polynomial efficiency) the operator R beg R tar . R beg is easy to compile; it's just 
a single multiply-controlled phase. Next, we will explain how to compile Rtar- 

2 Pcrhaps some symmetry of the physical situation being modeled by the Bayesian network x 
will suggest some x value that has non-zero probability. Alternatively, one can proceed as follows. 
For defmiteness, consider a Bayesian net x = (x.i,x 2 ,x 3 ) with 3 nodes. Suppose P{x3, X2, xi) = 
P(x3\x2, xi)P(x2\xi)P(xi) and the functions P^x^,^ Px 2 \x t an d Px x are known. Choose y\ £ Sx ± 
such that Px^TJi) 7^ 0. Then choose j/2 € Sx 2 such that Px 2 \x (2/2I2/1) 7^ 0. Finally, choose 2/3 G Sx 3 
such that Px^x^x^y^iVi) ^ 0. Set x = (2/1,2/2,2/3)- 

3 We will discuss in a future paper what to do if y/n(xo) is much larger than 0(l/\/Ns). 
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Figure 5: Definition of operator Vp for Szegedy operator W. 



FiglSl which is identical to Fig. 18 in Ref.jl], defines an operator Vp in terms 
of multiple (modified) phase estimation steps. The Vp of Ref. [3] depends on a param- 
eter (3 (inverse temperature) because the operator M in that paper depends on this 
parameter. Vp in the present paper does not depend on (3 (because the Bayesian net 
x doesn't) so we will drop the subscript from it henceforth, and refer to it simply 
as V . V does not depend on (3 but it still depends on the positive integers a and c. 
(In the language of Ref. [3], a=number of probe bits for each PE (Phase Estimation) 
step, and c=number of PE steps). Note that operator W is applied 2 a c times by V. 

Let |0 ac ) = \0)®( ac \ J = {0,±l,±2,...,±(iV 5 - 1)}, and J' = J - {0}. 
According to Lemma 2 of Ref. [3], for any e 2 > 0, if we adjust the integers a and c so 
that 

r 53 z = ° { 7i> ' (82) 

and 

c w log 2 (-^) , (83) 

V^2 



(recall S = 1 — |mi| is the distance between the two largest eigenvalue magnitudes of 
Mi), then V acts on the space Vb usy ® |0 ac ) as follows: 

|o-)(o-| v | Xi )(o-| 

v ~ \mu + h \i>m\ (84) 

where the \xj) are states of ac qubits such that (0 ac |xj) = 0. In Eq. (l8l|) and for the 
remainder of this section, the top row represents the ac ancilla qubits shown in FigEl 
and the bottom row represents the 2iY# qubits on which W operates. 



15 



Now define 



Q = (-l)l° ac >< oac l = l-2|O ac ><0° 



and 



Rtar = Q V . 



It follows that for any G Vb usy , 



R 



tar 



|0 ac ) 



|0 ac )(0 a 



1 _ 2V ] /v " 1 V 



|0 ac ) 



1 - 2 



|0 ac )(0 ac 

|-0O> <-0o| 

|0 ac ) |0 ac ) 



|0 ac ) 



+ 



|0 ac ) 

-Rtar I 



-2|^o)(^o|)|^) 



55) 



^6) 



57) 



(89) 
(90) 



Eg. ( 1901) is the essence of Corollary 2 in Ref. [3]. It means that R tar acting on Vbusy 
can be approximated by R tar acting on Vbusy <S> |0 ac ). Since we already know how to 
compile Rtar, we have accomplished our goal of compiling Rt ar , at least approximately. 

Next, we will try to estimate the error of our quantum Gibbs algorithm. Sup- 
pose 7f() is our estimate of 7r(). Note that for any x G S^, 



\tt(x) - 7f(x)| = |(a/7t(x) - ^/n(x))(y / n(x) + a/7t(x))| 

< 2\^{x) - ^{x)\ . 

Suppose e > is defined so that 



max | a/ tt(x) 

X 



'mx)\ < e . 



(91) 
(92) 



(93) 



Then, since we apply the Rb eg Rtar operator a total of L times, and each time we can 
incur an error of 



e « Lyfe . (94) 

If we define one step as one W application, then the total number of steps for the 
whole algorithm is 0(L2 a c) = 0(± log 2 (f )). Thus, our algorithm will yield a sample 

of the classical Bayesian net x with precision O(e), in C(^= log 2 (7)) steps. Achieving 
the same precision with a classical Gibbs sampling algorithm would require 0(j) 
steps. 
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The Szegedy operator W of this paper can also be used to do quantum sim- 
ulated annealing and Metropolis-Hastings if the marginals P(x* +1 |X| i | C ) can be cal- 
culated for each % from the transition matrix P(x* +1 |x*). (In the case of simulated 
annealing, P{x t+l \x t ) is different for each of the annealing schedule). 
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